version 18
set more off
quietly log
local logon = r(status)
if "`logon'" == "on" { 
	log close 
	}
log using bbw2025prq-analyze.log, text replace


/*		********************************************************************	*/
/*     	File Name:		bbw2025prq-analyze.do									*/
/*     	Date:   		January 14, 2025										*/
/*      Author: 		Frederick J. Boehmke									*/
/*      Purpose:		This file replicates the MLE of opinion with spatial  	*/
/*						decay reported in Table 1 and then interprets it by		*/
/*						creating Figures 2 and 3.								*/
/*      Input File:		bbw2025prq.dta											*/
/*      Output File:	bbw2025prq-analyze.log									*/
/*      				bbw2025prq-analyze-table1.txt							*/
/*      				bbw2025prq-analyze-figure2.png							*/
/*      				bbw2025prq-analyze-figure3.png							*/
/*		Requires:		estout.ado												*/
/*		********************************************************************	*/


	/*********************************************************************/
	/* Install required programs. Run the commented out lines if needed. */
	/*********************************************************************/

/*
net install estout, from(http://fmwww.bc.edu/RePEc/bocode/e)
*/

	/***********************************************************************/
	/* Define the two programs that implement the MLE for the ordered      */
	/* logit that estimates the spatial decay parameter. There are two     */
	/* versions fo the program: one for the ordered outcome variables with */
	/*three values and one for the outcome variable with five values.      */
	/***********************************************************************/

		/* Clear the two MLE programs so they can be reloaded when the file is re-run. */

capture program drop lf_decay3 lf_decay5

		/* MLE for outcomes with 3 values. */

program define lf_decay3

	args lnf Xb tau1 tau2 b_exposure lndelta

	tempvar pi1 pi2 exposure
	
		/* Re-generate the decayed distance from every casino */
		/* using with the current value of the decay parameter. */
		/* dist_miles1, dist_miles2, etc. contain the distances in miles. */
		/* Note that to ease estimation we paramaterize delta as its log so that */
		/* it can take on positive and negative values. Thus we use lndelta, */
		/* where delta = exp(lndelta) since delta = exp(ln(delta)). It is transformed back */
		/* to the original scale as delta when we create our final tables. */
	
	foreach var of varlist dist_miles* {
	
		quietly generat double `var'_decay = exp(-exp(`lndelta')*`var')
	
		}

		/* Add up the decayed distances over all casinos. */
		
	quietly egen double `exposure' = rowtotal(dist_miles*_decay)
	
		/* Calculate the probabilities for the first two outcome values. */
		
	quietly gen double `pi1' = invlogit( `tau1' - `Xb' - `b_exposure'*`exposure')
	quietly gen double `pi2' = invlogit( `tau2' - `Xb' - `b_exposure'*`exposure') - `pi1'

		/* Assign probabilities based on the observed outcome values. */
		
	quietly replace `lnf' = ln(`pi1') if $ML_y1==1
	quietly replace `lnf' = ln(`pi2') if $ML_y1==2
	quietly replace `lnf' = ln(1 - `pi1' - `pi2') if $ML_y1==3
	
		/* Drop the decayed distance variables since they will */
		/* need to be recalculated on the next iteration. */
	
	quietly drop dist_miles*_decay

end

		/* MLE for outcomes with 5 values. */

program define lf_decay5

	args lnf Xb tau1 tau2 tau3 tau4 b_exposure lndelta

	tempvar pi1 pi2 pi3 pi4 exposure
	
	foreach var of varlist dist_miles* {
	
		quietly generat `var'_decay = exp(-exp(`lndelta')*`var')
	
		}

	quietly egen `exposure' = rowtotal(dist_miles*_decay)
		
	quietly gen double `pi1' = invlogit( `tau1' - `Xb' - `b_exposure'*`exposure')
	quietly gen double `pi2' = invlogit( `tau2' - `Xb' - `b_exposure'*`exposure') - `pi1'
	quietly gen double `pi3' = invlogit( `tau3' - `Xb' - `b_exposure'*`exposure') - (`pi1' + `pi2')
	quietly gen double `pi4' = invlogit( `tau4' - `Xb' - `b_exposure'*`exposure') - (`pi1' + `pi2' + `pi3')

	quietly replace `lnf' = ln(`pi1') if $ML_y1==1
	quietly replace `lnf' = ln(`pi2') if $ML_y1==2
	quietly replace `lnf' = ln(`pi3') if $ML_y1==3
	quietly replace `lnf' = ln(`pi4') if $ML_y1==4
	quietly replace `lnf' = ln(1 - `pi1' - `pi2' - `pi3' - `pi4') if $ML_y1==5
	
	quietly drop dist_miles*_decay

end


	/*************************************************/
	/* Read in the data and summarize the variables. */
	/*************************************************/

	
use bbw2025prq, clear

		/* Only those in the post-election survey were asked these questions. */
		/* The values differ for 4 cases in the Common Content and the Iowa */
		/* content, but the relevant variables are all missing on the Iowa */
		/* version, so drop those extra 4 cases. */

tabulate tookpost tookpost_iow, missing
keep if tookpost_iow == 1

		/* Run the tabulations that are reported in Tables SI.1 through SI.3. */
		/* These use the full set of respondents asked these questions. */

tab1 ig_ben_trib_rev ig_ben_st_rev ig_cost_st_rev ig_opin_rev, missing

	
	/******************************************/
	/* Drop cases excluded from the analysis. */
	/******************************************/

	
drop if inlist(state, "AK", "HI", "DC")

		/* Remove cases with ideology = not sure. */

drop if ideo5==6 

		/* Create the summary statistics reported in Table SI.4. */

tabstat female age educ white natlec ideo5 raceresent fedtribe gaming gamnum /// 
	reli_bornagain reli_import reli_attend reli_protestant reli_catholic reli_other income, ///
	statistics(mean sd min max n) columns(statistics) format(%9.2f)
	
	
	/*********************************************************/
	/* Run the models reported in Table 1 in the main paper. */
	/*********************************************************/

		/*** Model 1: Benefits tribe. ***/
		
		/* Run an ordered logit without exposure to get starting values. */
	
ologit ig_ben_trib_rev female age educ white natlec ideo5 raceresent fedtribe gaming gamnum /// 
		  reli_bornagain reli_import reli_attend reli_protestant reli_catholic reli_other income
		  
		/* Save the coefficient vector to set up starting values. */
		  
matrix _b0 = e(b)
		
		/* Add starting values for the exposure coefficient and decay paramter. */
	
matrix _b0 = _b0 , 0.04 , `=ln(.02)'
		
		/* Specify and estimate the MLE model. */
	
ml model lf lf_decay3 (ig_ben_trib_rev = female age educ white natlec ideo5 raceresent fedtribe gaming gamnum /// 
		  reli_bornagain reli_import reli_attend reli_protestant reli_catholic reli_other income, noconstant) ///
		  /tau1 /tau2 /b_exposure /lndelta, vce(cluster stateno) ///
		  technique(nr)
		  
	ml init _b0, copy		  
	ml maximize, difficult 
	
		/* Store the results to access later. */
	
	estimates store ig_ben_tr
	
		/*** Model 2: Benefits state. ***/
	
ologit ig_ben_st_rev female age educ white natlec ideo5 raceresent fedtribe gaming gamnum /// 
		  reli_bornagain reli_import reli_attend reli_protestant reli_catholic reli_other income
		  
matrix _b0 = e(b)
		
matrix _b0 = _b0 , 0.04 , `=ln(.02)'
		
ml model lf lf_decay3 (ig_ben_st_rev = female age educ white natlec ideo5 raceresent fedtribe gaming gamnum /// 
		  reli_bornagain reli_import reli_attend reli_protestant reli_catholic reli_other income, noconstant) ///
		  /tau1 /tau2 /b_exposure /lndelta, vce(cluster stateno) ///
		  technique(nr)
		  
	ml init _b0, copy		  
	ml maximize, difficult 
	
	estimates store ig_ben_st
	
		/*** Model 3: Costs for state. ***/
	
ologit ig_cost_st_rev female age educ white natlec ideo5 raceresent fedtribe gaming gamnum /// 
		  reli_bornagain reli_import reli_attend reli_protestant reli_catholic reli_other income
		  
matrix _b0 = e(b)
		
matrix _b0 = _b0 , -0.04 , `=ln(.02)'
		
ml model lf lf_decay3 (ig_cost_st_rev = female age educ white natlec ideo5 raceresent fedtribe gaming gamnum /// 
		  reli_bornagain reli_import reli_attend reli_protestant reli_catholic reli_other income, noconstant) ///
		  /tau1 /tau2 /b_exposure /lndelta, vce(cluster stateno) ///
		  technique(nr)
		  
	ml init _b0, copy		  
	ml maximize, difficult 
	
	estimates store ig_cost_st
	
		/*** Model 4: Opinion on gaming. ***/
	
ologit ig_opin_rev female age educ white natlec ideo5 raceresent fedtribe gaming gamnum /// 
		  reli_bornagain reli_import reli_attend reli_protestant reli_catholic reli_other income
		  
	estimates store olog_ig_opin1

matrix _b0 = e(b)
		
matrix _b0 = _b0 , 0.02 , `=ln(.005)'
		
ml model lf lf_decay5 (ig_opin_rev = female age educ white natlec ideo5 raceresent fedtribe gaming gamnum /// 
		  reli_bornagain reli_import reli_attend reli_protestant reli_catholic reli_other income, noconstant) ///
		  /tau1 /tau2 /tau3 /tau4 /b_exposure /lndelta, vce(cluster stateno) ///
		  technique(nr)
		  
	ml init _b0, copy		  
	ml maximize, difficult 
	
	estimates store ig_opin1
	
		/* Test whether both coefficients are zero. */
		/* Zero for the decay makes the variable a constant, i.e., the number of gaming facilities. */
	
	testnl ([/]b_exposure=0) (exp([/]lndelta)=0)	

	lrtest olog_ig_opin1 ig_opin1, force
	
	
	/********************************************/
	/* Combine the results to generate Table 1. */
	/********************************************/

		/* Adjust variable labels to fit the table. */
		
label variable ig_opin 			"Opposition to Indian Gaming"
label variable ig_ben_trib 		"Gaming Does not Benefit Tribes"
label variable ig_ben_st 		"Gaming Does not Benefit State"
label variable ig_cost_st 		"Gaming Does not Cost State"
label variable ig_opin_rev 		"Support for Indian Gaming"
label variable ig_ben_trib_rev	"Gaming Benefits Tribes"
label variable ig_ben_st_rev	"Gaming Benefits State"
label variable ig_cost_st_rev	"Gaming Costs State"
label variable female 			"Female"
label variable age 				"Age"
label variable educ 			"Education"
label variable white 			"Anglo"
label variable natlec 			"National Economy"
label variable ideo5 			"Ideology"
label variable raceresent 		"Racial Resentment"
label variable fedtribe 		"Federal Tribes (#)"
label variable gaming 			"Indian Gaming"
label variable gamnum 			"Number of Gaming Tribes"
label variable reli_bornagain 	"Born Again"
label variable reli_import 		"Importance of Faith"
label variable reli_attend 		"Regularly Attend Church"
label variable reli_protestant 	"Protestant"
label variable reli_catholic 	"Catholic"
label variable reli_other 		"Other Religion"
label variable income			"Income"

		/* Output Table 1 to a text file. */
		/* Note that the decay parameter is transformed from the log scale */
		/* used during estimation back to the linear scale, i.e., delta = exp(lndelta). */
	
estout ig_ben_tr ig_ben_st ig_cost_st ig_opin1 using bbw2025prq-table1.txt, replace ///
	cells(b(star fmt(3)) se(par)) ///
	transform(lndelta exp(@) exp(@)) ///
	rename(tau1 cut1 tau2 cut2 tau3 cut3 tau4 cut4 b_exposure Exposure lndelta delta) ///
	order(eq1: cut*) ///
	stats(N ll, fmt(0 2)) ///
	modelwidth(8) ///
	starlevels(+ 0.10 * 0.05) legend ///
	label varwidth(20) ///
	collabels(, none) 
		

	/**********************************************************/
	/* Graph the distribution of distance to tribal gaming.   */
	/**********************************************************/
	
		/* Create a kernel density plot and histogram of closest. */

egen dist_miles_min = rowmin(dist_miles*)

kdensity dist_miles_min, scheme(s1mono) ///
	xlabel(#8) ///
	ylabel(#5, grid angle(0)) ///
	xtitle("Distance to Closest Casino") ///
	ytitle("Proportion of Respondents")	
	
hist dist_miles_min, scheme(s1mono) bin(32) fraction ///
	xlabel(#8) ///
	ylabel(#5, grid angle(0)) ///
	xtitle("Distance to Closest Casino") ///
	ytitle("Proportion of Respondents")
	
		/* How many respondents have a gaming facility within 50 or 100 miles? */
	
recode dist_miles_min (0/50 = 1)   (50/max=0), generate(within50)
recode dist_miles_min (0/100 = 1) (100/max=0), generate(within100)

tab1 within50 within100, missing

		
	/*******************************************************/
	/* Graph the estimated decay effects for each model.   */
	/*******************************************************/

		/* Load the decay parameter and exposure coefficient from each */
		/* model into scalars to include in the figures. */
	
foreach estimates in ig_opin1 ig_ben_tr ig_ben_st ig_cost_st {

	estimates restore `estimates'
	
	scalar delta_`estimates' = exp(/lndelta)
	scalar gamma_`estimates' = /b_exposure
	
	}

		/* Create Figure 2, which plots the estimated decay functions. */
	
twoway hist dist_miles_min, bin(32) fraction scheme(s1mono) yaxis(2) ///
	fcolor(gs12) lcolor(gs6) ///
 || function y = exp(-delta_ig_opin1*x), range(0 600) yaxis(1) ///
	lcolor(gs0) lpattern(shortdash) lwidth(medthick) ///
 ||	function y = exp(-delta_ig_ben_tr*x), range(0 600) ///
	lcolor(gs0) lpattern(solid) lwidth(medthick) ///
 ||	function y = exp(-delta_ig_ben_st*x), range(0 600) ///
	lcolor(gs0) lpattern(dash) lwidth(medthick) ///
 ||	function y = exp(-delta_ig_cost_st*x), range(0 600) ///
	lcolor(gs0) lpattern(longdash) lwidth(medthick) ///
    ylabel(none, axis(1)) yscale(range(0 1) axis(2)) xlabel(0(100)600) ///
    ylabel(0(0.25)1, grid axis(2)) ///
	xtitle(Distance of Respondent to Tribal Gaming Location, size(medlarge)) ///
	ytitle("Contribution of Tribal Gaming" "Location to Total Exposure", size(medlarge) axis(2)) ///
	ytitle("", axis(1)) ///
	legend(cols(1) label(2 Support for Gaming) label(3 Benefits Tribe) label(4 Benefits State)  ///
	  label(5 Costs State) ring(0) bplacement(0) order(2 3 4 5)) ///
	xsize(6) ysize(3)
	
	graph export bbw2025prq-analyze-figure2.png, replace width(3000)
	
		/* Create Figure 3, which multiplies the previous by the coefficient on exposure. */
	
twoway function y = gamma_ig_opin1*exp(-delta_ig_opin1*x), range(0 200) scheme(s1mono) ///
	lcolor(gs0) lpattern(shortdash) lwidth(medthick) ///
 ||	function y = gamma_ig_ben_tr*exp(-delta_ig_ben_tr*x), range(0 200) ///
	lcolor(gs0) lpattern(solid) lwidth(medthick) ///
 ||	function y = gamma_ig_ben_st*exp(-delta_ig_ben_st*x), range(0 200) ///
	lcolor(gs0) lpattern(dash) lwidth(medthick) ///
 ||	function y = gamma_ig_cost_st*exp(-delta_ig_cost_st*x), range(0 200) ///
	lcolor(gs0) lpattern(longdash) lwidth(medthick) ///
	ylabel(#5, grid) xlabel(0(50)200) ///
	xtitle(Distance of Respondent to Tribal Gaming Location, size(medlarge)) ///
	ytitle("Contribution of Tribal Gaming" "Location to Latent Opinion", size(medlarge)) ///
	legend(cols(1) label(1 Support for Gaming) label(2 Benefits Tribe)  label(3 Benefits State)  ///
	  label(4 Costs State) ring(0) bplacement(n)) ///
	xsize(6) ysize(3) 
	
	graph export bbw2025prq-analyze-figure3.png, replace width(3000)
	
clear
log close
exit, STATA  
